Introduction

This document shows how to calculate candidate Pillars and their constituent Indicators at Census Block Group, Tract, and Service Area/ Municipality level. Here, we load the data we pre-processed in the previous exercise

b <- read_sf("../data/boundary.geojson")
bg.l <- read_sf("../data/blockgroups_long.geojson")
bg.w <- read_sf("../data/blockgroups_wide.geojson")
tr.l <- read_sf("../data/tracts_long.geojson")
tr.w <- read_sf("../data/tracts_wide.geojson")
pl.w <- read_sf("../data/place_wide.geojson")
pl.l <- read_sf("../data/place_long.geojson")
lead <- read_sf("../data/parcels_lead_service_lines.geojson")
lead.c <- sf::st_centroid(lead)

calculate_volume <- function(hh_size = 4, gpcd = 50){
  vol = hh_size * gpcd * 30
  return(vol)  
}

calculate_bill <- function(volume = 6000){
  fixed_charges = 7.63 + 9.85 + 1.8
  vol_rate = 5.76 + 2.71 
  bill = fixed_charges+ (vol_rate*volume)/748.052
  return(bill)
}



gradeAffordability <-function(HBI,PPI){
 grade <- NA
 grade <- case_when(
    HBI < 7 & PPI < 20 ~ "Low Burden",
    HBI < 7 & PPI>=20 & PPI <35 ~ "Moderate-Low Burden",
    HBI < 7 & PPI >=35 ~ "Moderate-High Burden",
    HBI >= 7 & HBI < 10 & PPI >= 35 ~ "High Burden",
    HBI >= 7 & HBI < 10 & PPI >= 20 & PPI < 35 ~ "Moderate-High Burden",
    HBI >= 7 & HBI < 10 & PPI < 20 ~ "Moderate-Low Burden",
    HBI >= 10 & PPI >= 35 ~ "Very High Burden",
    HBI >= 10 & PPI >= 20 & PPI < 35 ~ "High Burden",
    HBI >= 10 & PPI < 20 ~ "Moderate-High Burden"
  )
  return(grade)
  
}


save.image("../cache/components_prep.RData")
load("../cache/components_prep.RData")

Pillar 1 | Access: Affordability and Service Quality

1.1 Affordability

1.1.1 HBI/PPI Matrix

We adopt the approach in the figure below, as elaborated in this report.

HBI/PPI Matrix

Overall (Census Place - level data)
#HBI Calculation

# Place
pl.w <- pl.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
pl.w <- pl.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
pl.w <- pl.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
pl.w <- pl.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
pl.w <- pl.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
pl.w <- pl.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
pl.w <- pl.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
pl.w <- pl.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))


pl.w$HBI_size_avg <- 100*pl.w$bill_size_avg/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_1 <- 100*pl.w$bill_size_1/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_2 <- 100*pl.w$bill_size_2/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_3 <- 100*pl.w$bill_size_3/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_4 <- 100*pl.w$bill_size_4/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_5 <- 100*pl.w$bill_size_5/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_6 <- 100*pl.w$bill_size_6/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_7 <- 100*pl.w$bill_size_7/(pl.w$hh_income_upper_limit_Q1E/12)


pl.w$PPI <- 100*(pl.w$inc_poverty_countE-pl.w$inc_poverty_count_gr2E)/pl.w$inc_poverty_countE
pl.w$aff_grade_avg <- gradeAffordability(pl.w$HBI_size_avg,pl.w$PPI)
pl.w$aff_grade_1 <- gradeAffordability(pl.w$HBI_size_1,pl.w$PPI)
pl.w$aff_grade_2 <- gradeAffordability(pl.w$HBI_size_2,pl.w$PPI)
pl.w$aff_grade_3 <- gradeAffordability(pl.w$HBI_size_3,pl.w$PPI)
pl.w$aff_grade_4 <- gradeAffordability(pl.w$HBI_size_4,pl.w$PPI)
pl.w$aff_grade_5 <- gradeAffordability(pl.w$HBI_size_5,pl.w$PPI)
pl.w$aff_grade_6 <- gradeAffordability(pl.w$HBI_size_6,pl.w$PPI)
pl.w$aff_grade_7 <- gradeAffordability(pl.w$HBI_size_7,pl.w$PPI)


pl.bill<- pl.w %>% st_drop_geometry() %>% select(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7) %>%
     pivot_longer(cols = c(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7),names_to="hh_size",names_prefix="bill_size_",values_to="bill")

pl.HBI<- pl.w %>% st_drop_geometry() %>% select(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7) %>%
     pivot_longer(cols = c(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7),names_to="hh_size",names_prefix="HBI_size_",values_to="HBI")

pl.aff<- pl.w %>% st_drop_geometry() %>% select(PPI,hh_income_upper_limit_Q1E,aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7) %>%
     pivot_longer(cols = c(aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7),names_to="hh_size",names_prefix="aff_grade_",values_to="aff_grade")

pl.aff <- pl.aff %>% left_join(pl.HBI,by="hh_size") %>% left_join(pl.bill,by="hh_size") %>% select(hh_size,bill,aff_grade,HBI,hh_income_upper_limit_Q1E,PPI) 


print(paste0("The Average Household Size for Naperville is ",
             round(pl.w$hh_size_avg_overallE,2)))

[1] “The Average Household Size for Naperville is 2.79”

print(paste0("At 50 GPCD, this corresponds to a monthly water volume of ",
             calculate_volume(hh_size=pl.w$hh_size_avg_overallE)))

[1] “At 50 GPCD, this corresponds to a monthly water volume of 4185”

print(paste0("This corresponds to an water and sewer bill of $",
             round(pl.w$bill_size_avg,2)))

[1] “This corresponds to an water and sewer bill of $66.67”

print(paste0("The Upper limit of the lowest quintile of annual household income for Naperville is $",
             pl.w$hh_income_upper_limit_Q1E,
             " or $",
             round(pl.w$hh_income_upper_limit_Q1E/12,2),
             " monthly"))

[1] “The Upper limit of the lowest quintile of annual household income for Naperville is $57958 or $4829.83 monthly”

print(paste0("The Poverty Prevalence Indicator (PPI) for Naperville overall is % ",
             round(pl.w$PPI,1)))

[1] “The Poverty Prevalence Indicator (PPI) for Naperville overall is % 9.4”

print(paste0("By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7"))

[1] “By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7”

datatable(pl.aff) %>%
  formatCurrency(c('bill','hh_income_upper_limit_Q1E'),'$') %>%
  formatRound(c('HBI','PPI'),2)

Below we visualize HBI, PPI, and the overall affordability grade as calcualted at the tract level.

# Tract
# Place
tr.w <- tr.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
tr.w <- tr.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
tr.w <- tr.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
tr.w <- tr.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
tr.w <- tr.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
tr.w <- tr.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
tr.w <- tr.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
tr.w <- tr.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))


tr.w$HBI_size_avg <- 100*tr.w$bill_size_avg/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_1 <- 100*tr.w$bill_size_1/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_2 <- 100*tr.w$bill_size_2/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_3 <- 100*tr.w$bill_size_3/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_4 <- 100*tr.w$bill_size_4/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_5 <- 100*tr.w$bill_size_5/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_6 <- 100*tr.w$bill_size_6/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_7 <- 100*tr.w$bill_size_7/(tr.w$hh_income_upper_limit_Q1E/12)


tr.w$PPI <- 100*(tr.w$inc_poverty_countE-tr.w$inc_poverty_count_gr2E)/tr.w$inc_poverty_countE
tr.w$aff_grade_avg <- gradeAffordability(tr.w$HBI_size_avg,tr.w$PPI)
tr.w$aff_grade_1 <- gradeAffordability(tr.w$HBI_size_1,tr.w$PPI)
tr.w$aff_grade_2 <- gradeAffordability(tr.w$HBI_size_2,tr.w$PPI)
tr.w$aff_grade_3 <- gradeAffordability(tr.w$HBI_size_3,tr.w$PPI)
tr.w$aff_grade_4 <- gradeAffordability(tr.w$HBI_size_4,tr.w$PPI)
tr.w$aff_grade_5 <- gradeAffordability(tr.w$HBI_size_5,tr.w$PPI)
tr.w$aff_grade_6 <- gradeAffordability(tr.w$HBI_size_6,tr.w$PPI)
tr.w$aff_grade_7 <- gradeAffordability(tr.w$HBI_size_7,tr.w$PPI)


hbi <- select(tr.w,
                  HBI_size_avg,
                  HBI_size_1,
                  HBI_size_2,
                  HBI_size_3,
                  HBI_size_4,
                  HBI_size_5,
                  HBI_size_6,
                  HBI_size_7,
                  )

ppi <- select(tr.w,
              PPI)

aff <- select(tr.w,
              aff_grade_avg,
              aff_grade_1,
              aff_grade_2,
              aff_grade_3,
              aff_grade_4,
              aff_grade_5,
              aff_grade_6,
              aff_grade_7)

m1 <- mapview(aff,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m2 <- mapview(hbi,burst=TRUE,col.regions=viridis::magma(7)) + mapview(b,alpha.regions=0,col.regions="green",stroke=TRUE,lwd=3,color="green",layer.name="Municipal Boundary")
m3 <- mapview(tr.w,zcol="hh_income_upper_limit_Q1E",col.regions=viridis::plasma(7),layer.name="Lower Quintile Income") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m4 <- mapview(ppi,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")


sync(m1,m2,m3,m4)

1.1.2, 1.1.3, 1.14, Cutoffs, Customer Assistance Programs (last 1 year), Conservation/ efficiency appliance incentive programs

Below, we simulate cutoff rates as a probabilistic function of HBI and PPI, and Customer Assistance participation rates as a probabilistic function of cutoff rates and white population share. We simulate participation in conservation/ efficiency appliance incentive programs as a function of homeownership rates and rate of population with at least a bachelors degree.

In reality, one would hope to use actually utility records of participation in these programs at the address level, or aggregated to census tracts.

tr.w$white_perc <- tr.w$pop_race_whiteE/tr.w$pop_race_countE

P1 <- select(tr.w, NAME, GEOID, hh_type_total_countE, aff_grade_avg,
             aff_grade_1,
             aff_grade_2,
             aff_grade_3,
             aff_grade_4,
             aff_grade_5,
             aff_grade_6,
             aff_grade_7,
             HBI_size_avg,
             HBI_size_1,
             HBI_size_2,
             HBI_size_3,
             HBI_size_4,
             HBI_size_5,
             HBI_size_6,
             HBI_size_7,PPI, hh_income_upper_limit_Q1E, hh_size_avg_overallE,white_perc, pop_owner_occupied_unitsE, pop_rental_unitsE, 
             pop_educ_attainment_countE, pop_educ_bachelorE, pop_educ_masterE, pop_educ_docE, pop_educ_profE) 

set.seed(20212002)

#Cutoffs and CAP
z <- 5  -0.03* P1$HBI_size_avg - 0.02*P1$PPI - 0.04 * P1$HBI_size_avg*P1$PPI
P1$Cutoff_Perc = 100*(1/(1+exp(z)))
P1$CAP_Perc <- P1$Cutoff_Perc * (1-0.3*(1-P1$white_perc))

#Conservation (ever)
P1$bachAbove <- 100*(P1$pop_educ_bachelorE + P1$pop_educ_masterE + P1$pop_educ_docE + P1$pop_educ_profE)/P1$pop_educ_attainment_countE
P1$own_rate <- 100*P1$pop_owner_occupied_unitsE/(P1$pop_rental_unitsE+P1$pop_owner_occupied_unitsE)

z <- 4 - 0.01 * P1$bachAbove - 0.02 * P1$own_rate
P1$Incentive_rate <- 100*(1/(1+exp(z)))

m5 <- mapview(P1,zcol="Cutoff_Perc", layer.name="Account cutoff (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m6 <- mapview(P1,zcol="CAP_Perc", layer.name="CAP participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m7 <- mapview(P1,zcol="Incentive_rate", layer.name="Efficiency Incentive participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")

sync(m5,m6,m7)

1.2 Water Quality (last 5 years)

We model all water quality events as basically random occurrences with probabilities between 0 and 10% and coverage of 50-100% at the tract level

1.2.1 Connections with Service Interruptions (not including nonpayment cutoffs or vacation-type shutoffs)

set.seed(32896)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$Interrup <- event*coverage

1.2.2 Connections with Boil water advisories

set.seed(328435396)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$boil_rate <- event*coverage

1.2.3 Health-based SDWA violations

This would be a Utility-Level measure only, and is available from EPA SDWIS

In this case, the API call for Naperville involves the PWSID IL0434670:

https://enviro.edap-cluster.com/cache/query/enviro3/efservice/VIOLATION/PWSID/IL0434670/IS_HEALTH_BASED_IND/Y/ROWS/0:100/JSON

In this case, there are no health-based violations in the past 10 years, so the value of the indicator is 0.

Below we export the Pillar 1 indicators.

P1_tract <- select(P1,-pop_owner_occupied_unitsE,-pop_educ_attainment_countE,-pop_rental_unitsE,-pop_educ_bachelorE,-pop_educ_masterE,-pop_educ_docE,-pop_educ_profE)
st_write(P1_tract,"../out/P1_tract.geojson",append=FALSE)
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.
## Deleting layer not supported by driver `GeoJSON'
## Deleting layer `P1_tract' failed
## Writing layer `P1_tract' to data source `../out/P1_tract.geojson' using driver `GeoJSON'
## Updating existing layer P1_tract
## Writing 45 features with 30 fields and geometry type Multi Polygon.
write.csv(st_drop_geometry(P1_tract),"../out/P1_tract.csv")


P1_overall <- P1_tract %>% st_drop_geometry() %>% summarise(across(c(HBI_size_avg:hh_income_upper_limit_Q1E,Incentive_rate:boil_rate,Cutoff_Perc:CAP_Perc),~weighted.mean(.,w=hh_type_total_countE))) %>% mutate(health_based_sdwa_violations=0)

write.csv(P1_overall,"../out/P1_overall.csv")

Pillar 2 | Economic and Workforce Development Contributions

Tract and Overall Level

  • Number Staff living in service area/tract
P2_tract <- tr.w %>% 
  select(NAME,GEOID,hh_type_total_countE) %>% 
  mutate(staff_count = round(runif(length(NAME),min=0,max=10),0),
         staff_per_1000hh = 1000*staff_count/hh_type_total_countE)

P2_overall <- P2_tract %>% st_drop_geometry %>% summarize(staff_count_in_service_area = sum(staff_count),
                                                          staff_per_1000hh = sum(staff_count)/sum(hh_type_total_countE))

P2_overall$staff_total<-300

Overall level Only

  • Diversity requirements in hiring

  • Diversity requirements in procurement

  • Local requirements for procurement/contracting

  • Participation in local education, training and workforce development programs

    • Binary
    • Trainees engaged with
set.seed(35455)
P2_overall <- P2_overall %>% 
  mutate(require_diverse_hiring = rbinom(1,1,0.5),
  require_diverse_procurement = rbinom(1,1,0.5),
  require_local_procurement = rbinom(1,1,0.5),
  regional_training_education_programs = round(runif(1, min=0,max=5)),
  region_training_education_participation = round(runif(1,min=0,max=100))
                       )

Below we export the Pillar 2 inidicators in a standardized format.

st_write(P2_tract,"../out/P2_tract.geojson",append=FALSE)
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.
## Deleting layer not supported by driver `GeoJSON'
## Deleting layer `P2_tract' failed
## Writing layer `P2_tract' to data source `../out/P2_tract.geojson' using driver `GeoJSON'
## Updating existing layer P2_tract
## Writing 45 features with 5 fields and geometry type Multi Polygon.
P2_tract <- st_drop_geometry(P2_tract)
write.csv(P2_tract,"../out/P2_tract.csv")
write.csv(P2_overall,"../out/P2_overall.csv")

Pillar 3: Resilience and Asset Management

Tract and Overall Level

  • Number Community water planning meetings (in each Tract)

  • Proportion structures with lead service lines (by tract)

  • Water main breaks/length pipe (by tract)

  • Sewer overflows (by tract)

  • Water quality complaints received/ addressed (by tract)

  • Lead service line replacements (proportion of lead service lines by tract)

  • Number leaks detected (by tract)

  • Expenditure incurred for detection measures

  • Infrastructure repair/rehabilitation/replacement expenditure (by tract)

  • Real losses reduced due to repairs/rehab/replacement (by tract)

P3 <- tr.w %>% 
  select(NAME,GEOID,hh_type_total_countE) 

P3.c <- st_join(P3,lead.c)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
P3.c <- P3.c %>% mutate(water_quality_complaints_received = rbinom(length(GEOID),1,0.08),
                    water_quality_complaints_resolved = water_quality_complaints_received*rbinom(length(GEOID),1,0.9),
                    ft_mains = runif(length(GEOID),55,164),
                    lead_line = as.numeric(SUBTYPE=="LEAD")) 

P3_tract <- P3.c %>% st_drop_geometry() %>%
  group_by(GEOID,NAME,hh_type_total_countE) %>% add_tally() %>% 
  summarize(conn = mean(n),
            water_quality_complaints_received = sum(water_quality_complaints_received),
            water_quality_complaints_resolved = sum(water_quality_complaints_resolved),
            ft_mains = sum(ft_mains),
            lead_lines = sum(lead_line,na.rm=TRUE)) %>%ungroup()
## `summarise()` has grouped output by 'GEOID', 'NAME'. You can override using the `.groups` argument.
P3_tract <- P3_tract %>% 
  mutate(meeting_count = round(runif(1, min=0,max=4)),
         meeting_count_1000hh = meeting_count/hh_type_total_countE,
         main_breaks = floor(runif(length(GEOID),min=0,max=3)),
         main_breaks_1000ft_pipe = main_breaks/(ft_mains/1000),
         prop_lead_lines = lead_lines/conn,
         leaks_detected = floor(runif(length(GEOID),min=main_breaks, max = 10)),
         prop_lead_lines_replaced = case_when(lead_lines > 0 ~ runif(1,min=0,max=1),
                                              lead_lines == 0 ~ NA_real_),
         sewer_overflows = floor(runif(length(GEOID), min=0, max=2)),
         leak_detection_expenditure = runif(length(GEOID),min=0,max=1000000),
         network_renewal_percent = runif(length(GEOID),min=0,max=20),
         network_renewal_expenditure =runif(length(GEOID),min=0,max=1000) * network_renewal_percent,
         real_losses_percent_delivered =  runif(length(GEOID),min=0,max=100),
         real_loss_reduction_percent = case_when(real_losses_percent_delivered > 5 ~ runif(length(GEOID),min=0,max=95),
                                                 real_losses_percent_delivered <= 5 ~ runif(length(GEOID),min=0,max=50))
         
         
         )

Overall Level Only

  • Utility Resilience Index (URI)

  • Infrastructure Leakage Index (by DMA -> summarize over Census Tract if possible)

P3_overall <- P3_tract %>% 
  summarize(hh = sum(hh_type_total_countE),
            conn = sum(conn),
            meeting_count = sum(meeting_count),
            main_breaks = sum(main_breaks),
            ft_mains = sum(ft_mains),
            lead_lines = sum(lead_lines),
            leaks_detected = sum(leaks_detected),
            lead_lines_replaced = mean(prop_lead_lines_replaced,na.rm=TRUE),
            sewer_overflows = sum(sewer_overflows)
           # network_renewal_percent2 = weighted.mean(x=network_renewal_percent,w=ft_mains,na.rm=TRUE),
            #real_losses_percent_delivered2 = weighted.mean(x=real_losses_percent_delivered,w=conn,na.rm=TRUE),
            #real_loss_reduction_percent2 = weighted.mean(x=real_losses_reduction_percent,w=conn,na.rm=TRUE)
            ) 
P3_overall <- P3_overall %>%  mutate(
            URI = floor(runif(1, min=0, max=100)),
            ILI = runif(1,min=1,max=4),
            prop_lead_lines = lead_lines/conn,
            breaks_per_1000ft = main_breaks/(ft_mains/1000),
            leaks_per_1000ft = leaks_detected/(ft_mains/1000),
            sewer_overflows_per_1000conn = sewer_overflows/conn
)

P3_overall$network_renewal_percent = weighted.mean(x=P3_tract$network_renewal_percent,w=P3_tract$ft_mains,na.rm=TRUE)
P3_overall$real_losses_percent_delivered = weighted.mean(x=P3_tract$real_losses_percent_delivered,w=P3_tract$conn,na.rm=TRUE)
P3_overall$real_loss_reduction_percent = weighted.mean(x=P3_tract$real_loss_reduction_percent,w=P3_tract$conn,na.rm=TRUE)

P3_tract <- left_join(P3,P3_tract,by=c("GEOID","NAME","hh_type_total_countE"))

Below we export the Pillar 3 indicators

st_write(P3_tract,"../out/P3_tract.geojson",append=FALSE)
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.
## Deleting layer not supported by driver `GeoJSON'
## Deleting layer `P3_tract' failed
## Writing layer `P3_tract' to data source `../out/P3_tract.geojson' using driver `GeoJSON'
## Updating existing layer P3_tract
## Writing 45 features with 21 fields and geometry type Multi Polygon.
P3_tract <- P3_tract %>% st_drop_geometry()
write.csv(P3_tract,"../out/P3_tract.csv")
write.csv(P3_overall,"../out/P3_overall.csv")

Overall Tract-level map

Below you can compare any two attributes side-by-side by selecting the appropriate layers in the layer control icon of each map in the upper left.

tracts <- left_join(P1_tract,P2_tract,by=c("GEOID","NAME"))
tracts <- left_join(tracts,P3_tract,by=c("GEOID","NAME"))
tracts <- select(tracts,-GEOID,-NAME)

t1 <- mapview(tracts,burst=TRUE,homebutton=FALSE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")

sync(t1,t1,ncol=1)

  1. ↩︎